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FOR  THERMAL  CONDUCTION  APPLICATIONS 
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Bradley  A.  Peavy 
ABSTRACT 


New  formulas  for  calculating  thermal  response  factors  for  mult iple -layer 
construction  have  been  developed  by  a rigorous  derivation.  A comparison 
was  made  of  the  time  for  computation  between  the  presently  used  matrix 
algebra  method  and  the  method  given  in  this  paper.  Results  were  obtained 
using  the  new  method  in  one-fiftieth  to  one-half  of  the  computational 
time  necessary  to  obtain  solutions  from  the  matrix  algebra  method. 

Comparisons  with  another  analytical  method  were  performed  to  verify  the 
accuracy  of  the  response-factor  technique. 
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DETERMINATION  AND  VERIFICATION  OF  THERMAL  RESPONSE  FACTORS 
FOR  THERMAL  CONDUCTION  APPLICATIONS 

by 

Bradley  A.  Peavy 

1 . INTRODUCTION 

The  analytical  solutions  needed  to  describe  temperature  or  heat  flow  for 
steady  periodic  or  transient  conduction  heat  transfer  in  multi-layer  walls, 
ceilings  and  floors  are  quite  complicated  and  are  not  available  for  non- 
linear conditions  such  as  thermal  radiation  at  surfaces  and  time-dependent 
changes  in  surface  film  coefficients.  It  has  therefore  been  expedient 
to  employ  "approximation"  methods  by  which  non-linear  conditions  may  be 
satisfied.  One  such  method,  which  is  the  result  of  an  analytical  formu- 
lation, is  termed  the  "response-factor  method,"  and  solves  for  temperatures 
or  heat  flows  of  multi-layer  constructions  based  on  the  past  temperature 
history . 

A response-factor  method  defined  by  Kusuda  [1]*  uses  overlapping  triangular 
pulses  to  compute  response  factors  for  a particular  construction.  These 
response  factors  are  then  used  to  determine  temperatures  and  heat  flows 
in  response  to  previously  occurring  events.  In  order  to  handle  multi-layer 
constructions  for  solution  of  response  factors  by  computers,  a matrix 
algebra  was  developed  for  an  arbitrary  number  of  layers.  There  are  some 
inherent  difficulties  in  using  this  approach,  in  that  the  calculation 
of  response  factors  can  be  quite  time  consuming  because  certain  functions 


* See  references  cited  at  end  of  text 
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are  not  well  suited  for  efficient  calculation.  A portion  of  this  paper 
will  be  devoted  to  the  development  of  equations  which  allow  for  increased 
computational  efficiency.  The  algorithms  presented  are  sufficient  for  the 
determination  of  up  to  a seven -layer  composite,  which  is  probably  more 
than  sufficient  for  possible  building  constructions. 

Unfortunately,  the  accuracy  of  the  temperatures  and  heat  flows  values  cal- 
culated by  response-factor  methods  has  received  very  little  verification. 
Although  it  is  well  known  that  the  thermal  properties  of  building  materials 
are  not  as  well  defined  as  they  should  be,  this  is  no  excuse  for  not  know- 
ing the  error  incurred  by  the  use  of  an  approximate  method  for  whatever 
thermal  properties  may  be  assigned  to  a particular  program  for  solution. 

A portion  of  this  paper  will  be  devoted  to  a comparison  between  numerical 
results  from  an  analytical  solution  and  results  from  calculations  using 
the  response-factor  technique. 


2.  ANALYSIS 

For  one-dimensional  heat  flow  in  an  individual  layer  of  one  or  more  parallel 
layers,  the  partial  differential  equation  is  given  by 
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where  vm  is  the  temperature  potential  above  a datum  plane,  x is  a dimen- 
sion along  which  heat  is  flowing,  is  the  thermal  diffusivity  of  the  layer 
material,  and  t is  the  time.  For  continuity  of  temperature  and  heat  flow 
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between  layers,  perfect  contact  is  assumed,  i.e., 
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where  K is  the  thermal  conductivity 
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the  Laplace  transform  to  (1)  gives 
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for  which  a solution  is: 
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where  b - is  the  distance  from  x=0  (surface  of  first  layer)  to  the 
m-1 

nearest  face  of  the  m-th  layer.  Applying  the  continuity  conditions  (2) 
gives  the  following  relations 
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and  is  the  thickness  of  the  (m-l)-th  layer.  This  process  is  continued 

until  the  constants  and  (pertaining  to  the  layer  with  face  exposed 
at  x=0)  are  found  in  terms  of  and  (pertaining  to  the  n-th  layer  with 
face  exposed  at  x=b  ) . At  x=0,  the  heat  flux  is  proportional  to  the 
temperature  difference  between  the  fluid  (air,  gas  or  liquid)  and  the 
surface,  and  is  represented  by  the  relationship 
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(4) 


where  is  the  surface  film  resistance  and  f(t)  is  fluid  temperature 
as  a function  of  time.  Similarly,  the  boundary  condition  at  x=b^  is 
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where  R0  is  the  surface  film  resistance  and  g(t)  is  a fluid  temperature 
as  a function  of  time.  When  either  R^  or  R£  is  zero,  the  time  temperature 
function  represents  temperatures  at  the  respective  surfaces.  The  resulting 
expressions  for  the  transform  of  the  temperatures  in  layer  1 and  layer  n 
are  given  by 
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(The  above  summations  are  over  m=l,2,...  2 , n being  the  number  of  layers.) 
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and  J , L and  A.  are  defined  in  Table  1. 
m m i,m 
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[The  transforms  of  the  heat  flux  at  x=0  and  at  x=bn  are  found  by  differen- 
tiating (6)  and  (7)  with  respect  to  x and  multiplying  by  minus  one  and  the 
rrespective  thermal  conductivity: 
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The  inversion  of  (8)  and  (9)  is  performed  by  evaluating  the  residues  at 
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the  poles  of  the  denominator  W=0,  where  p=  -8  or  v/p  = i8 , which  gives 
the  relationship 
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The  residues  at  the  poles  p=  -3  are 
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Of  particular  concern  is  the  evaluation  of  the  first  root  of  (10) . This 
can  be  done  expeditiously  by  expanding  the  sines  and  cosines  in  their 
series  and  considering  only  the  first  two  terms,  in  order  to  obtain  an 
initial  estimate  of  the  first  root 
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Consider  the  triangular  pulse  function  of  Kusuda  [1],  where 
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which  when  substituted  for  f(p)  and  g(p)  in  (8)  and  (9)  gives  double 
poles  at  p=0.  Following  are  limits  of  the  necessary  functions  of  (8) 
and  (9)  for  evaluating  the  residues  at  the  double  poles. 
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For  the  first  term  in  (8),  the  residue  for  0<t£6  is 


K1  |~tBl  A1B2  A2B1 

6/^  L A1  A2 


(16) 


9 


for  the  last  term  in  (9) , the  residue  is 
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and  for  the  first  term  in  (9)  and  the  last  term  in  (8) , the  residue  is 
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and  R is  the  total  thermal  resistance  of  the  n-layers  plus  R^  and  R2 
The  response  factors  X^,  and  evaluated  at  t=6,  become 
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For  Y 2 and  Z ^ evaluated  at  t=26. 
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For  larger  values  of  j , the  response  factors  decrease  with  increase  in 


j by  a common  ratio  — i.e.,  for  N sufficiently  large, 
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for  all  j>N,  where  3^  is  the  first  root  of  (10).  Thus  for  j greater  than 
N,  the  response  factors  can  be  computed  by  the  relationship 


The  temperatures  f(t)  and  g(t)  may  be  constructed  from  a series  of  triangular 
pulse  functions  (15),  overlapping  in  time  by  an  amount  6.  By  the  principle 
of  superposition,  the  heat  flux  at  time  T is  determined  from  the  sum  of  the 
products  of  the  response  factors  and  the  temperatures  from  the  present  time 
to  preceding  time  intervals  of  duration  6.  The  heat  flow  at  time  t=T  and 
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at  the  x=0  face  is  then 
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and  the  heat  flux  at  x=b  is 
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where  V is  the  temperature  potential  in  relation  to  a fixed  datum  plane 
temperature,  which  is  referenced  to  time,  t=T , T-l,  t-2,  etc. 


(24) 


For  certain  types  of  constructions,  the  number  of  required  response 

factors  can  become  quite  large  to  give  a reasonable  accuracy  to  the  heat 

fluxes  of  (23)  and  (24).  Conduction  transfer  functions  can  reduce  this 

number  considerably  and  are  defined  from  the  relationship,  F - CR  F 

1 ,T  1 , 

of  (23)  and  similarly  from  *,24). The  conduction  transfer  functions  are 
then  defined  as  follows: 
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X-  = X--CR  X-  , 
J J J-l 


Using  these  functions, 


Yj = YrCR  h-i  aj  ■ arCR  aj-i 


the  heat  flux  at  t=T  and  x=0  becomes 


(25) 


Fl,x  ■ CR  • - YjVn,x-j*l> 

J=1 


From  (21)  and  (25),  it  can  be  seen  that  the  conduction  transfer  functions 
are  numerically  very  small  numbers  as  j approaches  N,  and  the  number  of 
functions  needed  for  computation  purposes  is  considerably  reduced.  However, 
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This  must  be 


initially  it  is  necessary  to  know  the  value  of  F^  T_^ 
determined  from  several  iterations  over  the  past  temperature  history 
of  and  Vn . 

3.  RESPONSE  FACTORS 

Using  the  algorithms  of  the  previous  section,  the  response  factors  of 
equations  (19),  (20),  and  (21)  are  calculated  from  the  computer  program 
found  in  Appendix  A.  To  show  time  savings  in  response-factor  calculations, 
identical  multi-layer  constructions  were  input  into  both  the  program  of 
Kusuda  [1]  and  the  program  of  Appendix  A.  Six  multi-layer  constructions 
are  shown  in  Table  2,  and  the  computation  times  for  the  two  programs  are 
shown  in  Table  3.  As  can  be  seen,  considerable  time  saving  is  available 
from  the  program  of  Appendix  A. 

In  the  example  of  Table  2,  there  is  an  evident  lack  of  enclosed  air  spaces 
in  the  building  construction.  Kusuda  [2]  assumed  only  a purely  thermal 
resistance  effect  of  air  spaces  during  dynamic  temperature  changes.  It 
is  the  supposition  of  the  author  that  both  the  air  space  thickness  and 
the  heat  capacity  of  the  air  should  be  considered  for  heat  transfer  cal- 
culations. From  literature  values,  the  thermal  diffusivity  of  dry  air 
varies  from  0.639  ft^/h  at  0 °F  to  .977  ft^/h  at  120  °F,  and  these  values 
are  reduced  somewhat  by  the  presence  of  water  vapor.  For  the  program  of 
Appendix  A,  the  thermal  diffusivity  of  air  in  air  spaces  is  assumed  to 

r\ 

be  0.75  ft^/h,  the  air  space  thickness  is  assumed  to  be  one  inch  if  not 
specified  by  the  input  card,  and  the  thermal  resistance  is  as  defined  for 
steady-state  air-space  values. 
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Table  2.  Multi-Layer  Constructions 


Layer 

Thermal 

Specific 

Layer 

Th  ickness 

Conductivity 

Density 

Heat 

Description 

ft 

Btu/h  ft  F 

lb/ft3 

Btu/lb  F 

Case  1 Two  Layers 


1 1-IN  MINERAL  FIBERBOARD 

0.833 

.035 

23.0 

.140 

2 4-IN  LIGHT-WEIGHT  CONCRETE 

.3333 

.100 

40.0 

.200 

Case  2 Three  Layers 

1 5/ 8-IN  PLASTER  BOARD 

.0521 

.094 

50.0 

.260 

2 3-1/2  IN  BATT  INSULATION 

.2917 

.026 

2.0 

.220 

3 3/4-IN  WOOD  SIDING 

.0625 

.800 

36.0 

.280 

Case  3 Four  Layers 

1 1/2-IN  PLASTER  BOARD 

.0417 

.094 

50.0 

.260 

2 3/ 4tIN  POLYSTYRENE  INSULATION 

.0625  • 

.094 

50.0 

.260 

3 4-IN  COMMON  BRICK 

.3333 

.417 

120.0 

.190 

4 4-IN  FACE  BRICK 

.3333 

.750 

130.0 

.190 

Case  4 Five  Layers 

1 1/2-IN  PLASTER  BOARD 

.0417 

.094 

50.0 

.260 

2 1-IN  BATT  INSULATION 

.0833 

.026 

2.0 

.220 

3 4-IN  L.W.  CONCRETE 

.3333 

.100 

40.0 

.200 

4 3/4-IN  BOARD  INSULATION 

.0625 

.017 

2.2 

.290 

3 1-IN  STUCCO 

.0833 

.400 

116.0 

.200 

Case  5 Six  Layers 

1 3/4-IN  ACOUSTIC  TILE 

.0625 

.035 

30.0 

.200 

2 4-IN  HEAVY  WEIGHT  CONCRETE 

.3333 

1.000 

140.0 

.200 

3 2-IN  ROOF  INSULATION 

.1667 

.025 

5.7 

.200 

4 2-IN  H.W.  CONCRETE 

.1667 

1.000 

140.0 

.200 

5 3/ 8-IN  FELT 

.0313 

.110 

70.0 

.400 

6 1/2-IN  SLAG 

Case  6 Seven  Layers 

.0417 

.830 

55.0 

.400 

1 

1/2-IN  PLASTER  BOARD 

.0417 

.094 

50.0 

.260 

2 

1-IN  BATT  INSULATION 

.0833 

.026 

2.0 

.220 

3 

1/2-IN  PLYWOOD 

.0417 

.067 

34.0 

.290 

4 

3-IN  H.W.  CONCRETE 

.2500 

1.000 

140.0 

.200 

5 

1/2-IN  PLYWOOD 

.0417 

.067 

34.0 

.290 

6 

3/4-IN  BOARD  INSULATION 

.0625 

.017 

2.2 

.29  0 

7 

1-IN  STUCCO 

.0833 

.400 

116.0 

.200 
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Table  3.  Execution  Time  for  Computation  of  Response  Factors,  Seconds 


Kusuda  [ 1 ] 


Appendix  A Time  Savings 


Two  layers 

2.037 

.187 

1.850 

Three  layers 

8.885 

.174 

8.711 

Four  layers 

1.9.02 

.378 

1.5  24 

Five  layers 

3.018 

.625 

2.39  3 

Six  layers 

1.609 

.846 

.763 

Seven  layers 

9.786 

1.681 

8.105 

Table  4. 

Heat  Flux  at 

Inside  Surface  of 

Wood  Frame 

O 

Construction  with  Air  Space,  Btu/h  ft 


Time  (hr) 

No  heat  capacity 

3 1/2"  air  space 

5 1/2"  air  space 

1 

1.742 

1.021 

.954 

2 

1 . 262 

1.242 

1.203 

3 

1.359 

1.342 

1.319 

4 

1.548 

1.525 

1 .494 

5 

1.707 

1.691 

1.666 

6 

.964 

1.109 

1 . 233 

7 

- 1.966 

- 1.598 

- 1.154 

8 

- 6.113 

- 5.684 

- 5.024 

9 

-10.410 

- 9.924 

- 9.214 

10 

-14.248 

-13.809 

-13.144 

11 

-17.427 

-17.065 

-16.501 

12 

-19.720 

-19.463 

-19.042 

13 

-20.908 

-20.783 

-20.544 

14 

-20.827 

-20.849 

-20.822 

15 

-19.497 

-19.660 

-19.846 

16 

-17.038 

-17.334 

-17.718 

17 

-13.592 

-13.996 

-14.552 

18 

- 9.416 

- 9.908 

-10.599 

19 

- 4.960 

- 5.459 

- 6.201 

20 

- 1.785 

- 2.088 

- 2.639 

21 

- .746 

- .870 

- 1.128 

22 

- .214 

- .281 

- .414 

23 

+ .239 

+ .185 

+ .093 

24 

.680 

.628 

.54  7 
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The  introduction  of  these  changes  does  alter  the  response  factors  and 
heat  flux  amplitudes  in  response  to  external  temperature  variations. 
Table  4 gives  heat  flux  values  at  the  inside  surface  for  a wood-frame 
construction,  with  no  heat  capacity  [2],  and  31/2  and  5 1/2  inch  air 
spaces.  The  external  temperature  variation  is  taken  from  Figure  1. 


4.  VERIFICATION  OF  RESPONSE  FACTORS 

Response  factors  as  defined  in  this  paper  are  the  result  of  an  analytical 
formulation  which  employs  a past  time  history  of  overlapping  pulses  for 
the  temperatures  at  or  adjacent  to  surfaces  of  composite  solids  to  deter- 
mine the  heat  flux  or  temperature  at  present  time.  A linear  variation 
in  temperature  over  the  time  period,  (usually  one  hour)  is  assumed.  The 
ability  of  response  factors  to  give  accurate  values  for  heat  flux  has 
been  questioned.  For  this  reason,  it  is  appropriate  to  compare  the  numeri- 
cal results  from  the  response -factor  calculation  with  those  from  another 
analytical  calculation. 

Analytical  solutions  can  be  found  from  equations  (6)  and  (7),  where  the 
temperature -t ime  functions  can  be  defined  by  the  trigonometric  series, 

An  p + BnWn 

f(t)  = E (An  cos  Wnt  + Bn  sin  Wnt ) ; f(p)  = £ r (27) 

p + W 
r n 


and 


g(t)  = CQ  + Z (Cn  cos  WRt  + Dn  sin  WRt ) 


g(p)  = — +2 
P 


Cp  + DW 
nr  n n 

2 2 
p + W 


(28) 
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TIME  OF  DAY  - HOURS 

Figure  1.  Temperature  variations  for  g(t)  and  f(t)  used  in  verification  of  response  factors. 


Considering  only  the  steady  periodic  condition,  the  residues  at  the  poles 
at  p=0  and  p=+iwn  can  be  found  for  the  surface  temperatures.  Heat  flow 
at  the  surfaces  can  then  be  obtained  from  (4)  and  (5). 

The  six  examples  of  Table  2 were  used  in  determination  of  heat  flux  for 
both  the  response  factor  and  analytical  solution.  The  temperature  varia- 
tion, g(t)  shown  in  Figure  1 was  used  for  the  test  cases,  with  f(t)=0 
and  R-^  = 0.68  and  R2  = .333  hft^F/Btu.  The  coefficients  CQ , CR  and  Dn 
were  determined  from  96  points  where  there  was  linear  interpolation 
between  the  hourly  points.  This  was  necessary  to  assure  a nearly  linear 
variation  in  temperature  between  the  hourly  points  for  the  function  g(t). 
When  coefficents  were  determined  for  24  and  48  points  the  agreement  between 
the  heat  flux  for  the  analytical  and  response  factor  was  not  as  good. 

Heat  flux  was  computed  both  from  the  response  factors  [(23)  and  (24)]  and 
the  conduction  transfer  function  (26).  For  all  cases  examined  the  heat 
flux  at  the  inside  surface  was  less  than  one  percent  different  from  the 
values  computed  for  the  analytical  solution.  At  the  outside  surface,  the 
agreement  was  not  as  good,  particularly  around  5 and  19  hours,  where 
there  are  sudden  changes  in  the  slope  of  the  temperature  curve  (as 
shown  on  Figure  1).  The  detail  figures  show  that  the  function  g(t)  used 
in  the  analytical  solution  varies  considerably  from  the  linear  form 
required  by  the  response  - factor  method  in  the  time  period  4 to  5 and  19 
to  20.  It  is  expected  that  if  more  points  were  used,  there  would  be 
a better  agreement  at  the  outside  surface. 
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Table  5 

. Comparison  of  Heat 

Flux  Values 

Calculated 

from 

Response 

-Factor  and 

Analyt ical 

Methods  Btu/h 

ftz 

Case 

2. 

Case 

3. 

Outside 

Surface 

Inside  Surface 

Outside 

Sur  face 

In  s i d e 

Sur  face 

Resp . 

Anal . 

Resp . 

Anal  . 

Resp . 

Anal  . 

Resp . 

Anal  . 

1.128 

1.103 

.358 

.358 

39.450 

39.437 

-5.468 

-5.468 

.487 

.503 

.442 

.442 

34.607 

34.609 

-4.867 

-4.866 

1.177 

1.152 

.479 

.478 

32.679 

32.668 

-4.295 

-4.295 

1.268 

1.195 

.542 

.542 

30.798 

30.768 

-3.763 

-3.763 

.652 

.304 

.604 

.604 

27.171 

27.036 

-3.271 

-3.270 

- 9.247 

- 9.338 

.430 

.434 

- 2.856 

- 2.920 

-2.816 

-2.816 

-15.479 

-15.366 

- .489 

- .486 

- 38.739 

- 38.726 

-2.408 

-2.408 

-16.531 

-16.412 

-1  .933 

-1.933 

- 65.508 

- 65.481 

-2.110 

-2.109 

-16.683 

-16.552 

-3.456 

-3.456 

- 84.931 

- 84.961 

-2.012 

-2.011 

-15.951 

-15.807 

-4.861 

-4.862 

- 96.946 

- 96.894 

-2.170 

-2.169 

-14.261 

-14.108 

-6.047 

-6.048 

-101 .138 

-101.076 

-2.590 

-2.589 

-11.537 

-11.379 

-6.930 

-6.932 

- 96.975 

- 96.904 

-3.238 

-3.237 

- 7.701 

- 7.565 

-7.432 

-7.434 

- 83.896 

- 83.828 

-4.058 

-4.057 

- 3.379 

- 3.286 

-7.488 

-7.490 

- 63.289 

- 63.233 

-4.978 

-4.978 

.709 

.783 

-7.091 

-7.093 

- 37.984 

- 37.933 

-5.916 

-5.916 

5.175 

5.183 

-6.284 

-6.284 

- 8.151 

- 8.124 

-6.786 

-6.787 

8.621 

8.604 

-5.110 

-5.111 

21.992 

22.007 

-7.509 

-7.510 

12.274 

12.076 

-3.661 

-3.662 

53.149 

53.093 

-8.015 

-8.016 

11.223 

10.820 

-2.065 

-2.064 

71 .912 

71.763 

-8.254 

-8.255 

1.296 

1.332 

- .811 

- .806 

58.711 

58.703 

-8.195 

-8.197 

1.226 

1.242 

- .333 

- .332 

52.360 

52.354 

-7.857 

-7.858 

1.366 

1.371 

- .113 

- .113 

48.294 

48.290 

-7.332 

-7.332 

1.520 

1.517 

.055 

.055 

45.399 

■ 45.394 

-6.722 

-6.722 

1 .676 

1.648 

.214 

.214 

43.178 

43.164 

-6.901 

-6.091 
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Comparison  of  calculated  heat  flux  values  from  the  analytical  and  response- 
factor  methods  is  shown  in  Table  5 for  Cases  2 and  3. 

5.  CONCLUSIONS 

Formulas  for  calculating  thermal  response  factors  for  plane  multi-layer 
constructions  have  been  developed  as  given  by  equations  (19),  (20),  and 
(21).  A computer  program  to  obtain  response  factors  based  on  these 
formulas  is  found  in  the  Appendix.  A comparison  was  made  of  the  time 
for  computation  of  response  factors  between  the  matrix  algebra  method 
of  Kusuda  [1]  and  the  method  given  in  this  paper.  A considerable  saving 
in  computation  time  is  realized  by  the  methods  of  this  paper. 

Presently,  the  response- factor  calculations  for  constructions  with  air 
spaces  assume  only  a purely  thermal  resistance  effect  of  air  spaces 
during  dynamic  temperature  variations  [2].  The  heat  transfer  across  an 
air  space  involves  the  nature  of  the  bounding  surfaces,  the  intervening 
air,  orientation  of  the  space  and  the  direction  of  heat  flow:  and  hence 
the  three  modes  of  heat  transfer — radiation,  convection,  and  conduction  — 
influence  heat  flow  in  an  air  space.  It  would  be  impractical  to  simu- 
late the  three  modes  of  heat  transfer  for  constantly  changing  air  space 
surface  temperatures,  but  it  is  felt  that  a reasonable  approximation 
should  include  both  the  heat  capacity  and  air  space  thickness  for  dynamic 
heat  transfer  calculations.  These  were  included  in  the  computer  program 
of  the  Appendix.  The  introduction  of  these  changes  gives  different 
values  for  the  response  factors  and  resulting  heat-flux  amplitudes. 
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Response  factors  are  analytical  formulations  from  the  partial  differential 
equations  for  heat  conduction.  When  properly  applied,  the  response- 
factor  method  gives  correct  values  for  temperature  and  heat  flow  for 
conduction  heat  transfer  problems. 
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7.  CONVERSION  FACTORS  TO  METRIC  (S.I.)  UNITS 


Physical  Quantity 

To  Convert 
From 

To 

Multiply  By 

Length 

ft 

m 

3.0480  E-l 

Area 

ft2 

2 

m 

9.2903  E-2 

Temperature 

F 

C 

( F—  3 2 ) / 1 .8 

Density 

lbm/ft3 

kg/m3 

1.6018  E+l 

Thermal  Conductivity 

Btu/h  ftF 

W/mK 

1.7296  E+0 

Thermal  Resistance 

h ft2F/Btu 

m2K/W 

1.7623  E-l 

Thermal  Diffusivity 

f t2/h 

m2/ s 

2.9900  E-3 

Heat  Flux 

Btu/h  ft2 

W/m2 

3.1525  E+O 

Specific  Heat 

Btu/lbmF 

j/kgK 

4.1840  E+3 

Time 

h 

s 

3.6000  E+3 

21 


APPENDIX  A 


Appendix  A gives  a Fortran  listing  of  computations  for  the  response 
factors  and  includes  input  of  data  identical  to  that  of  Kusuda  [1]. 
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C PROGRAM  PCS  COMPUTING  RESPONSE  FACTORS  OF  MULTI-LATER  BUILDING 

C CONSTRUCTIONS  (ONE  TO  SEVEN  LAYERS > FOR  THERMAL  CONDUCTION 

C APPLICATIONS 

C INPUT 

C DEL-  TIMB  INCREMENT  ( EQU  15) 

C JA-  NUMBER  OF  LAYERS  IN  CONSTRUCTION  (MAY  INCLUDE  SURFACE  FILMS) 

C THERMOPHYSICAL  PROPERTIES  OF  EACH  LAYER  (ONE  CARD  PER  LAYER- 1 OF 7 . 0 ) 

C SURFACE  FILM  RESISTANCES  MAT  BE  INCLUDED  AS  LAYERS-  IN  'RE1  ONLY 
C IF  INCLUDED  FOR  BOTH  SIDES  - 'JA1  SHOULD  NOT  EXCEED  9 

C ORDERING  OF  LAYERS  FROM  INSIDE  TO  OUTSIDE 
C L - THICKNESS  OF  EACH  LAYER 

C X - THERMAL  CONDUCTIVITY  OF  LAYER 

C D - DENSITY  OF  LAYER 

C SP-  SPECIFIC  HEAT  eF  LAYER 

C RE-  THERMAL  RESISTANCE  (FOR  AIR  SPACES  OR  SURFACE  FILMS  ONLY) 

C DESCRIPTION  OF  EACH  LATER  (ONE  CARD  PER  LAYER-6A6  ) 

C RM-  CONTAINS  DESCRIPTIONS  FOB  EACH  LAYER  (HOLLERITH) 

C BLANK  CARD  INDICATES  NO  MORE  CONSTRUCTIONS  TO  BE  INPUTTED 

C 

C 

C SUBROUTINE  ROOTS  - COMPUTES  ROOTS  OF  EQU  10 

C 

C SUBROUTINE  ABC  - COMPUTES  (1-1)  TERMS  FOR  BQUS  16.17,16 

C (1-0)  SUMMATION  TERMS  FOR  EQUS  10,11 

C 

c 

C OUTPUT 

C THERMOPHYSICAL  PROPERTIES  OF  EACH  LAYER  AND  DESCRIPTIONS 

C U N ).K(  N),D(  N),SP(  N ).RE(  N >.  RM(  I,N  ),  1-1,6  (EXCLUDING  SURFACE  FILMS) 

C STATEMENT  - SURFACE  FILMS  INCLUDED  OR  EXCLUDED 

C X - PRINTOUT  OF  TERMS  USED  IN  COMPUTING  EQUS  16.17,18 

C BB, BC.BD.BE.BF.BG  - TERMS  OF  BQUS  16.17,16 

C Y - RCOTS  OF  EQU  10 

C CB-  COMMON  RATIO  EQU  22 

C RESPONSE  FACTORS 

C XX(  N ),YY(  N ).ZZ(  N),  N-l.IB  IB  DEFINED  BY  EQU  22 

C 

PARAMETER  F-60.G-100 

DIMENSION  S(  64  ),B(  7 ),C(  6 ),B(  6 ),L(6  ).K(  6 ).D(  8),SP(  6 ) , RE  ( 8),RM(  6.6), 
1AU  7 ).X(  F ).Y(  F ).W(  F ).T<  F ).XX(  G ).YT(G  ),ZZ(  G ) 

COMMON/CBA/S.C 
REAL  X.L 

100  PORMAT  ( 1017) 

102  FORMAT  ( 10F7.0) 

105  FORMAT  ( 6A6 ) 

201  FORMAT  ( 10F12.6  ) 

202  FORMAT  (3F16.6) 

203  FORMAT  ( 1 1 0.F1 1 . 4, FI  0. 3. F10 . 1 . FI 0 .3. F8. 2. 2I.6A6  ) 

205  FORMAT  ( 26H  BUILDING  CONSTRUCTION  NO. 14  ) 

206  FORMAT  ( 33B  ROOTS  OF  CHARACTERISTIC  EQUATION) 

207  FORMAT  ( 24E  RESPONSE  FACTORS  X.Y.Z) 

206  FORMAT!  1 HI  ) 

209  FORMAT!  1 H ) 

210  FORMAT  ( 56 E RESPONSE  FACTORS  ARE  CALCULATED  FROM  SURFACE  TO  SURFAC 
AB  ) 

211  FORMAT  ( 56 B RESPONSE  FACTORS  INCLUDE  SURFACE  FILM  RESISTANCES  - R1 
A-F6.3.8H  AND  R2-F6.3) 

212  FORMAT  (14B  COMMON  RATI0-F9.6) 

READ  < 5,102  ) DEL 

IR-0 
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97  WRITE  <6. 206) 

99  ULEAD  ( 5,100)  JA 
J “JA 
81  -.0 
82“. 0 

IF  ( J.BQ.O  ) 06  TO  200 

READ  (5.102)  L(  1 ).K(  1 ).D(  1 ),  £P(  1 ).  BE(  1 ) 

I A “2 

IF  < D( 1 ).GT. .001  ) GO  TO  98 
81 -BE( 1 ) 

J-J-l 
I A - 1 

96  DO  101  N-IA.J 

101  BEAD  (5.102)  U N ).  K(  N ).  D(  N ) , SP(  N ) . BE(  N ) 
I B- IB “1 

-READ  (5,105)  ( BM(  1.1  ),  1-1,6  ) 

103  DO  10*  N-IA.J 

10*  BEAD  (5.105)  ( EM(  N,  I ). 1-1,6) 

IF  (D(J).GT.  .001  ) 00  TO  106 
B2-BECJ) 

J-J-l 

106  IF  ( J.GT.7.0B. J.EQ.O  ) GO  TO  200 
DO  1 1 0 N-l ,J 

IF  ( D(  N ) ) 109,107.109 

107  AL(  N )-  • 75 

IF  ( L(  N ).LT.  .001  ) L(N)-. 083333 
E(  N )-L(  N )/BE(  N ) 

GO  TO  110 

109  AL(  N)-K(  N)/(D(N)*SP(  N)> 

110  CONTINUE 

AA-.O 

19  DO  20  N-l. 7 

B(  N )-. 0 

IF  (N.LE.J)  AA-AA-U  N )/X(  N ) 

20  IF  (N.LE.J)  B(  N )-U  N )/SQBT(  AL(  N ) ) 

AD-1  . 

IF  ( J.EQ.l  ) GO  TO  1 
I-J-l 

DO  21  N-l  , I 

E(  N )-I(  N-l  )*SQET(  AL(  N )/AL(  N-l  ) )/*(  N ) 

C<  N )-(  1 . -E(  N ) )/(  1 . -E(  N ) ) 

21  AD-2.*AD/(  1 . *E(  N ) ) 

GO  TO  ( 1 .2.3.*. 5.6.7  ).J 

7 S<  32  )-B<  7 )-B(  6 )*B(  5 )-B(  * )*B(  3 )-B(  2 )-B(  1 ) 
8(  31  )-B<  7 )-B(  6 )-B(  5 )-B(  * ) -B(  3 )*B(  2 )-B(  1 ) 
S<  3 0 )-B(  7>-B<  6 )-B(  5)-B(  ♦ )-B(  3 )*B(  2 )*B(  1 ) 
8(  29  )-B(  7 >-B<  6 )-B(  5 )-B<  * )-B<  3 )-B(  2 )-B(  1 ) 
8(2  8 )-B(  7 )-B(  6 )-B(  5 )-B(  * )-B(  3 )-B(  2 )-B(  1 ) 
S(  27  )-B(  7 )-B<  6 ) «B(  5 )-B(  * )-B(  3 )-B(  2 )-B(  1 ) 
S(  26  >-B<  7 )-B(  6 )-B(  5 )»B(  * >-B(  3 )-B(  2 )-B(  1 ) 
S(  25  )-B<  7 )-B(  6 )-B(  5 ) -B(  * )-B(  3 )-B(  2 )-B(  l ) 
S(  2*  ) -B(  7 )-B(  6 )-B(  5 )*B(  ♦ ) -B(  3 )-B(  2 )*B(  1 ) 
S(  23  )-B(  7 )-B(  6 )-B(  5)-B(  * )-B(  3 l-B(2  )-B(  1 ) 
S(  22  )-B<  7 )— B(  6 )-B(  5 >-B<  ♦ )-B(  3 )-B(  2 )-B(  1 ) 
S(  21  )-B(  7 )-B<  6 )-B(  5 )-B(  * )-B(  3 )-B(  2 )*B(  1 ) 
S(20)-B(  7 )-B(  6 )-B(  5 )-B(  ♦ )-B(  3 )-B(  2 )*B(  1 ) 
S(  19)-B<  7 )-B(  6 )-B(  5 )-B(  * )-B(  3 )-B(  2 )-B(  1 ) 
S<  1 8 )-B(  7 )-B<  6 )-B(  5 )-B(  * )-B(  3 )-B(  2 )-B(  1 ) 
S(  17  )-B(  7 )-B(  6 )-B(  5 )-B(  ♦ )-B(  3 )-B(  2 )-B(  l ) 

6 S(  16)-B(  7 )*B(  6 )-B(  5 )-B(  * )-B(  3 )-B(  2 )-B<  1 ) 
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S(  15  l-BC  7 l-BC  6 l-BC  5 l-BC  * l-BC  3 l-BC  2 l-BC  1 > 

S ( 1*1  -BC  7 I -B(  6 )-B(  5 1 *B(  * 1*B(  3 l-BC  2 l-BC  1 ) 

6(  1 3 1 -B<  7 ) -BC  6 l-BC  5 l-BC  * l-BC  3 l-BC  2 1 -BC  1 1 
6<  12  l-BC  7 l-BC  6 l-BC  5 )-B(  ♦ 1-B(  3 l-BC  2 l-BC  1 1 

S(  1 1 l-BC  7 1 -BC  6 l-BC  S l-BC  * l-BC  3 l-BC  2 l-BC  1 1 

SC  10  1 -BC  7 1 -BC  6 l-BC  5 l-BC  * l-BC  3 l-BC  2 l-BC  1 1 

6C  9 1-B( 7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 1 
5 SC  8 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 1 

SC  7 l-BC  7 l-BC  6 l-BC  5 l-BC  ♦ l-BC  3 1 -BC  2 1 -BC  1 1 

SC  6 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 1 

SC  5 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 1 

4 SC  4 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 1 

SC  3 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 I 

3 SC  2 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 I 

2 SC  1 l-BC  7 l-BC  6 l-BC  5 l-BC  4 l-BC  3 l-BC  2 l-BC  1 I 

GO  TO  9 
1 SC  1 l-BC 1 1 
9 DO  8 H-1,32 
I-N-32 

8 SC  I 1-SC  N 1-2. *BC  1 1 
AB-KC  1 1/SQBTC  ALC  1 I 1 
AC-EC  J 1/SQBTC  ALC  J 1 1 
T1-E1*AB 
V2-B2*AC 
I-l 

CALL  ABCC CA.X.J.f  1 
CA-T1*V2 
CB-T2 -VI 
CC-V2-T1 

BA-XC  2 1-XC  6 1-CB*XC  1 l-CC*XC  5 1 
BB -AB*C  XC  1 1-XC  5 1 1/BA 
BC-AC*C  XC  1 1-XC  5 1 1/BA 
BD-AC*AD/BA 

BH-C  XC  4 1-XC  8 1 1/6.  -CA4C  XC  2 1-XC  6 1 1-C  CB*XC  3 1-CC*XC  7 11/2. 

BI-C  XC  3 1-XC  7 1 1/2.  -V2*C  XC  2 1-XC  6 1 1 

BJ-C  XC  3 1-XC  71  1/2.  -VI *C  XC  2 1-XC  61  1 

BE-AB*C  BA*BI-BB*C  XC  1 1-XC  5 111/C  DEL*BA*BA  1 

BP- AC*C  BA*BJ-BH*C  XC  1 1-XC  5 1 1 1/C  DEL*BA*BA  1 

BG--AC*AD«BH/  C DEL*BA*BA  1 

XC  1 1-SQBTC  BA/BH  1 

CALL  BeOTSC AB.AC.AE.V1.V2.DEL.T.X.T.W. J.M 1 

IXC  1 1-BB-BE 

ZZC  1 l-BC -BP 

TTC  1 1-BD-BG 

XXC  2 1— BE 

ZZC  2 1— BP 

TTC  2 1--BC 

DO  50  N-l.M 

XXC  1 1-XXC  1 1-XC  N 1 
ZZC  1 1-ZZC  1 1— TC  N 1 
TTC 1 l-TTCl 1-TC W 1 
CA -EXPC -DEL*TC  N 1**2 1-2. 

XXC  2 1-XXC  2 1-XC  N 1*CA 
ZZC  2 l-ZZC  2 1-TC  N 1*CA 
5 0 TTC  2 1-TTC  2 1-TC  K1*CA 
CB-EXPC  -DEL*TC  1 1**2  1 
DO  56  1-3. G 
XXC  1 1-.0 
TTC  I l-.O 
ZZC  I l-.O 
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IB-1 

CC-I-3 

CA-CC*Y( 1 >**2*DEL 
IF  (CA.GT.20.)  GO  TO  59 
DO  55  N-l.M 

CA-DEL*T(  N)4*2 
CB-CC*CA 

IF  (CB.GT.40.)  GO  TO  56 
CD-BXP(  -CA  ) 

CD -EXF(  -CB>*(  l.-CD  )**2 
IKI  )-ZX(  I >-X<  N )* CD 
ZZ(  I >-ZZ<  I >-V(  N >*CD 

55  Tt(  I )-YT(  I )-T(N)*CD 

52  CA -ABS(  XXC  I )/XX(  1-1  )-CB  )*ABS(TT(  I >/YY(  1-1  >-CB  ) 

CA  -CA  *ABS(  ZZ(  I )/ZZ(  1-1  )-CB  ) 

56  IF  ( CA.LT.. 00003  ) GO  TO  59 
59  VBITE  (6.205)  IB 

DO  20*  I-l.J 

20*  VBITE  ( 6.203)  I . L(  I ).  X(  I ).  D(  I ).  SP(  I ).  BE(  I ),  ( BN(  I.N  ),N- 1 . 6 ) 
VBITE  (6.209) 

IF  (J.EQ.JA)  VBITE  (6.210) 

IF  (J.NE.JA)  VBITE  (6.211)  B1.B2 
VBITE  (6.206) 

VBITE  (6.201)  (TOO.N-1.5) 

VBITE  (6.209) 

VBITE  (6.207) 

VBITE  (6,202)  ( XX<  N ),  TTC  N ),  ZZ(  N ),  N-l  , I B ) 

VBITE  (6,212)  Cfi 
GO  TO  97 
200  STOP 
END 


SUBROUTINE  ABC(X.Z.J.I) 

DIMENSION  2X  1 ).T(  6*  ),V(64.4  ) 

COMMON/CBA/SC  64  ),C(  6 ) 

X-l 

IF  ( I.EQ.l  ) GO  TO  20 

M-l 

IF  (J.GT.l)  M-244J/4 
DO  6 N-l.M 
DO  6 1-1.2 
L-N 

IF  ( I.BQ.2  ) L-N-32 
B-X«S(  L ) 

A-EIN(B) 

B-COSI  B ) 

▼<  L , 1 ) -A 
L.2  )-B 

Y(L.3)-B*S(  L) 

8 ▼(  L.4  )-A«S(  L ) 

12  DO  14  N-l. 64 
1*  T(N)-V(N.X) 

15  T-.O 

v-.o 

10  GO  TO  (1. 2, 3. 4. 5. 6. 7). J 

7 A-C(  1 >*<  C(  2 )*T(  24  )-C(  3 )«T(  26  )♦«  4 )*(  T(  27)-C(  2 )*C(  3)*T(  32  ) ) )*T(  17  ) 
A— C(  5 )•(  A-C<  2 )*<  Ci  3 )*T(  29  )*C(  4 >*T<  30  ) )*C(  3 )«C<  4 )«T(  31  ) >♦«  2 )*T(  20) 
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B-C<  * HKT<  18  )*C(  2 >*<  Ci  1 >*T(  23  )*<X  3 >*T<  28)  >-Ot  1 >•«  3 )«TC  25)  ) 
A-A-B-C(  3 )*<T<  IS)  -CC  1 >*C(  2 )*7<  22  ) )*C<  1 )#T(  21  ) 

Y-Y-CI  6)*A 

A«C<  2 )*(C(  3 )*C<  4 )*TC  60  )®C(  3 )*CX  5 )♦«  61  )*C<4  )*C(  5)*T(  62  ) ) 

B-C(  2 )«T(  52  )*C(  3 )*T<  51  )*C<  4 >*T<  50>-C<  5 )•<  T(  49  )♦«  3 )*Ct  4 )*TC  63  ) > 
B-T(  53 )-C(  1 )*(  A*B  )-C<  3 )*C  C(  4 )*T<  57>-C<  5 )*T(  58  ))*C(4  )*«  5 )*T<  59  ) 
A-C(  2 )*(  C(  3 )*T<  54  ) *C(  4 )»T(  55  ) *C<  5 )*(  T(  56  )♦«  3 )*CX  4 )«T(  64  ))  ) 
WW-C<  6 )*{  A-B  ) 

6 A-C<  4 )*<  T(  9 )*C(  1 )*(  £X  2 )*T(  1 4 ) 3 )*T<  15)0*02  )*T(  11  )♦«  1 )*T<  12  ) 

Y-Y-C<  5 )*(  A*CC  3 >*<  T(  10  )-C(  2 >*<  C<  1 )*T(  13  )♦«  4 >*«  16)  ) ) ) 

A-K44  )*C(  1 )*<  C(  4 )*T(  41  )-C(  3 )*T(  42  )*C(  2 )*T(  43  ) )♦«  3 )•&  4 )*T<  47  ) 
W-W-C(  5)»<  A-C(  2 )*<C(  3)*T{  45  )*CX  4 )*(  K 46)*0  1 )»C(  3 )#7<4  8 > ) ) ) 

5 Y-Y-CI4  )*(  C<  1 )4T(  5 )*C(  2 )*T<  6)*0  3)*(T(7  )*C<  1 )*C<  2 >*T<  8)  ) ) 

»-W-C(  4 )*(  T(  37  )*C(  1 )*<  C(  2 )*T(  38)»«3  )*T(  39  ) )♦«  2 )*«  3 )*T(40  ) ) 

4 Y -Y  *C<  3 )*(  C<  1 )*T(  3 >*C<  2 )*T(  4 ) ) 
l-W*C(3  )*(  T(  35  )*C<  1 )*C<  2 )*T<  36)  ) 

3 Y-Y*C(  1 )*C(  2 )*7<  2 ) 

B-W*C(  2)«T(  34  ) 

2 Y-Y-7C  1 ) 

»-W-C(  1 )*K33  ) 

ZC  X >-Y 
Z(  Z *4  ) -» 

16  1-1*1 

IF  ( I.BQ.l  ) 06  76  22 
IF  ( K.LE.4  ) GO  76  12 
KETOBN 
1 Z(  I ) -T( 1 > 

IX  K*4  )-.0 
06  76  16 

20  D6  21  N-1,64 

21  T<N)-1. 

G6  76  15 

22  DG  23  N-1,64 

23  T<  N )-T(  N )*S(  N ) 

IF  ( K.LE.4  ) 06  76  15 

BETGSN 

BND 


S0BB6CTINE  B66TSC  AB , AC, AD, V 1 , V2 , BEL, B ,Z. P , fi, J , X ) 
DIMENSION  B(  1 ),Z4  1 ).P(  1 ).B(  1 ),Y(  8 ) 

1-0 

V-V14V2 

G-V2-Y1 

H-V2-V1 

F-2. 

U-.02 

M-l 

K-0 

A-.005 

CALL  ABC(A.Y.J.I) 

D "Y(  1 )*Y(  5 )*A*{  G*Y(  2 )*H*Y(  6 )-A*V*C  YC  1 )-Y<  5 ) ) ) 

X-ZC  1 ) 

E°.l 5*X 

4 CALL  ABC(X.Y,J,I) 

A -YC  1 )-Y(5)*X*<G*Y(2)*B*Y(6)-X*V*<Y(l  )-Y(  5 ) ) ) 

T-(  l.-Y*I*I)»Y(3)-(l.-V*I#X  )*Y(  7>-2.*X*V*(  Y(  1 )-Y<  5 ) ) 
IF  <ABS(  A).LT.l  ,E-6)  GO  76  7 
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IF  ( K.GT.2  ) 03  TO  6 

C INCREMENTING  TO  OBTAIN  FIRST  ZERO  - AMD  APPROXIMATE  FOR  OTHERS 
IF  ( A*D  ) 5.5.8 

5 X-X-E 
E "E/2. 

IF  ( M.BQ.l  ) GO  TO  4 
K-K*l 
GO  TO  4 
8 X*X ♦£ 

GO  TO  4 

C MEW TON "RAPE SON  METBOD  FOR  SOLVING  BQU  10 

6 Q-T  «G*<  J(  2 >-X*T(  4 ) I^BwC  T<  6 )-X*T(  8 ) ) 

X-X-A/Q 

K-K*l 

C NOT  TO  EXCEED  8-TB  ITERATION 
IF  C E.GT.8  ) GO  TO  7 
GO  TO  4 
T B(  M >-X 

A *EXF( -X*X*DEL  ) 

Q-T*G*<  T(  2 l-X*T(4  J )♦  B*(  T<  6 )-X»T(  8 ) ) 

E -2  . *A/<  X*  X*DEL*0  ) 

Z(  M )-E*<  Y(  2 )-T(  6 )-X*V2*<  T<  1 )-T<  5 ) ) HAB 
P( M )-E*AC*AD 

R(  M )«E«AC*(  Y(  2 )»T(  6 )-X*Vl*(  T<  1 )-T<  5 ) ) ) 

IF  ( X«X*DEL.GT.30. ) RETURN 
IF  (M.EQ.60)  RETURN 

X-0 

IF  ( M.GT.l  ) GO  TO  1 
IF  ( X.GT..S)  GO  TO  1 
F-3.5 
U-.002 

1 x-x»u 

CALL  ABC( X.T.  J.  I ) 

D-T(  1 )*Y<  5 )*X*(  G*T(  2 )*B*T(  6 )-X«Vw(TC  1 >-Y<  5 )>  ) 

IF  ( M.EQ.l  ) E-X/4. 

IF  < M.GT.l  I B-(  X-B(  M-l  ) )/F 

M"M  *1 

X"*X*E 

GO  TO  4 

END 
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